adams_moulton_coefficients(1, 1)
adams_moulton_coefficients(1, 2)
adams_moulton_coefficients(2, 3)
adams_moulton_coefficients(3, 4)
adams_moulton_coefficients(4, 5)

function adams_moulton_coefficients(s, p)
    % 定义符号变量 b0, b1, ..., b_{s}
    syms b [1 s+1]  % 创建 b1, b2, ..., b{s+1} (MATLAB下标从1开始，b1对应b0)

    % 设置阶条件方程
    eqns = sym(zeros(1, p));  % 需要 p 个方程
    for m = 1:p
        j = 0:s;  % j = 0, 1, ..., s
        if m == 1
            eqns(1) = sum(b) == 1;  % C1 条件: b0 + b1 + ... + b_s = 1
        else
            % 高阶条件: sum(j.^(m-1) * b_j) = [s^m - (s-1)^m] / m
            eqns(m) = sum(j.^(m-1) .* b) == (s^m - (s-1)^m)/m;
        end
    end

    % 求解方程组
    sol = solve(eqns, b);

    % 提取并显示系数（转换为分数形式）
    fprintf('\nAdams-Moulton 系数 (s=%d, p=%d):\n', s, p);

    % Handle different solution formats
    if s == 1 && p == 1
        % Special case when s=1 and p=1
        fprintf('b0 = %s\n', char(rats(double(sol.b1))));
    elseif s+1 == 1
        % When only one coefficient exists
        fprintf('b0 = %s\n', char(rats(double(sol))));
    else
        % General case with multiple coefficients
        for k = 1:s+1
            bk = sol.(sprintf('b%d', k));
            fprintf('b%d = %s\n', k-1, char(rats(double(bk))));
        end
    end
end